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Predicting  the  liquid  film  dynamics  inside  the  injector  cup  of  gas-centered  swirl  coax¬ 
ial  fuel  injectors  requires  a  general  two-phase  approach  that  is  appropriate  for  all  liquid 
volume  fractions,  high  Weber  number,  and  complex  geometries.  The  rapid  exchange  of 
momentum  at  the  highly  convoluted  interface  requires  tight  numerical  coupling  between 
the  gas  and  liquid  phases.  An  Eulerian  two-phase  model  is  implemented  to  represent  the 
liquid  and  gas  interactions  in  the  injector  as  well  as  the  atomization  processes  at  the  rough 
interface.  The  model,  originally  proposed  by  Vallet  et  al,  assumes  that  in  the  limit  of 
infinite  Reynolds  and  Weber  number,  features  of  the  atomization  process  acting  at  large 
length  scales  are  separable  from  small  scale  mechanisms.  A  transport  equation  for  the 
liquid  volume  fraction  represents  the  dispersion  of  the  liquid  into  the  gas  via  a  traditional 
turbulent  diffusion  hypothesis.  A  model  for  the  growth  of  mean  interfacial  surface  area  is 
then  used  to  characterize  the  growth  of  instability  at  the  interface,  allowing  a  characteriza¬ 
tion  of  Sauter  mean  diameter.  The  model  shows  promise  as  a  computationally  inexpensive 
tool  for  characterizing  spray  quality  in  regions  where  optical  experimental  data  are  difficult 
to  obtain  and  two-phase  direct  numerical  simulation  methods  are  too  demanding.  Two- 
dimensional  results  of  the  model  are  compared  to  photographs  of  the  liquid  film  within  the 
injector  cup.  The  comparisons  show  good  agreement  between  the  predicted  film  profile 
and  experimental  measurements,  although  in  some  cases  underpredicting  the  film  length. 
The  discrepancy  between  the  experiment  and  model  results  suggest  the  need  to  extend  the 
inter-phase  coupling  to  a  form  more  suitable  for  anisotropic,  swirling  flow. 


Nomenclature 


0  Farve  averaged  quantity 

0  volume  averaged  quantity 

0'  fluctuating  quantity 

k  turbulent  kinetic  energy 

e  dissipation  rate 

Ht  turbulent  viscosity 

Ui  mean  velocity  components 

p  mean  pressure 

p  mean  density 

Y  mean  liquid  mass  fraction 

Y  mean  liquid  volume  fraction 

S  mean  interfacial  surface  area  density 

a  surface  tension  coefficient 
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Introduction 


Fuel  injection  elements  in  modern  rocket  propulsion  systems  operate  at  extremely  high  Reynolds  and 
Weber  number,  producing  a  wide  range  of  relevant  length  scales.  Rendering  these  disparate  length  scales 
with  direct  numerical  simulation  techniques  is  extremely  expensive,  since  the  mesh  resolution  would  need  to 
be  fine  enough  to  resolve  the  smallest  droplets  yet  have  an  extent  that  spans  the  entire  nozzle.  Experimental 
techniques,  though  able  to  provide  macroscopic  insights,  have  difficulty  quantifying  the  small  scale  features 
of  the  spray  formation  due  to  the  primary  atomization  zone’s  topological  complexity  and  optically  dense 
nature. 

Gas-centered  swirl-coaxial  injectors  (GCSC)  are  particularly  challenging  to  characterize  due  to  the  fact 
that  atomization  occurs  within  a  small,  visually  obscured  cup  and  is  accompanied  by  a  fully  three-dimensional 
flow.  The  injector  functions  by  injecting  high  velocity  gaseous  oxidizer  axially  through  the  center  of  the 
element,  while  liquid  propellant  is  injected  tangentially  along  the  element  wall  to  produce  a  swirling  liquid 
film.  A  spray  is  formed  by  a  combination  of  instability  growth  along  the  surface  of  the  film  and  shear 
induced  by  the  high  velocity  co-flowing  gas.  Although  the  injector  resembles  other  injection  schemes  such  as 
pressure-swirl  and  coaxial  injectors,  the  GCSC  differs  in  that  the  atomization  occurs  within  the  injector  cup 
before  the  liquid  enters  the  combustion  chamber,  without  forming  a  conical  sheet  under  typical  operating 
conditions.  To  effectively  predict  atomization  quality,  it  is  essential  both  to  predict  the  uniformity  of  the 
spray  and  to  accurately  predict  film  length.  For  designs  in  which  the  injector  cup  is  significantly  longer 
than  the  film  length,  overheating  of  the  injector  faceplate  can  occur.  For  cup  lengths  shorter  than  the  film, 
incomplete  atomization  may  occur. ^ 

In  an  extension  of  their  2006  paper  in  which  they  outline  general  mechanisms  by  which  atomization  takes 
place,  Lightfoot  et  al.  published  a  report  in  which  they  determine  the  relevant  processes  at  play  within  the 
GCSC.  They  suggest,  based  on  a  inviscid  linear  instability  analysis,^  that  liquid  instability  growth  cannot 
be  a  dominant  mechanism  based  on  the  prediction  that  few  waves  would  reach  a  height  capable  of  causing 
atomization,  and  that  they  would  only  be  capable  of  atomizing  approximately  2.5%  of  the  total  film  volume. 
Under  similar  arguments,  they  also  suggest  that  liquid-phase  turbulence  cannot  be  a  dominating  factor,  and 
that  the  dominant  mechanism  taking  place  is  the  impact  of  gas-phase  turbulent  structures  with  the  liquid 
interface. 

Based  on  this  assumption,  they  propose  a  model  by  which  cylindrical  gas  eddies  “scoop”  out  a  disturbance 
in  the  liquid.  Although  the  majority  of  the  assumptions  made  are  only  to  give  a  rough  characterization  of 
entrainment  rate,  they  propose  a  useful  criterion  for  gas  turbulence  driven  atomization  to  occur:  the  kinetic 
energy  of  the  turbulent  eddy  must  exceed  the  surface  tension  energy  of  the  interface. 

Pg^eddy  —  (^) 

^eddy 

They  use  this  equation  to  develop  a  criteria  to  judge  what  fraction  of  eddies  are  large  enough  to  cause 
stripping,  and  use  turbulence  data  from  fully  developed  pipe  flow  to  estimate  an  eddy  size  distribution. 
This  relation  is  also  useful  for  verifying  the  assumption  that  gas-phase  turbulence  is  dominant.  Under 
the  assumption  that  the  dominant  eddies  lie  within  the  range  of  the  larger,  energy  containing  scale,  the 
characteristic  length  and  velocity  scales  are^ 


hi 

L  =  —  (2) 

C/  =  fc5  (3) 

The  assumption  that  deddy  and  v^ddy  are  on  the  order  of  L  and  U  respectively,  gives  the  following  relation 
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Lightfoot  et  al.  later  performed  experiments  upon  a  variety  of  GCSC  geometries  to  obtain  measurements 
of  film  length  as  a  function  of  relevant  non-dimensional  groups,  such  as  momentum  flux  ratio  and  mixture 
ratio.  Figure  1  demonstrates  the  range  of  geometry  considered  and  a  sample  correlation. 

Lightfoot’s  work  provides  a  quantitative  characterization  of  the  film  within  the  injector  cup.  Downstream 
measurements  of  atomization  quality  have  been  performed  over  a  wide  variety  of  operating  conditions  over 


2  of  13 


American  Institute  of  Aeronautics  and  Astronautics 


20 

15 

^10 

5 

0 


i  odhdtu  ♦onpntn  4-ouputu 
■  odpdtd  Bonpdtd  aouhdtd 


0.5  1  1.5  2  2.5  3  3.5 


ml/in2 


Name 

To  (mm) 

■c  (mm) 

r,  (mm) 

s  (mm) 

ODHUTD 

7.620 

1.321 

3.429 

2.870 

ODPDTD 

7.620 

1.321 

5.461 

0.838 

ODHNTN 

7.620 

1.651 

4.445 

1.524 

ODPDTN 

7.620 

1.651 

5.461 

0.508 

ODHUTU 

7.620 

1.981 

3.429 

2.210 

ODHDTU 

7.620 

1.981 

5.461 

0.178 

ONPDTD 

9.525 

1.321 

5.461 

2.743 

ONHNTD 

9.525 

1.321 

6.350 

1.854 

ONPDTN 

9.525 

1.651 

5.461 

2.413 

ONPNTN 

9.525 

1.651 

6.350 

1.524 

ONPUTN 

9.525 

1.651 

7.468 

0.406 

ONPNTU 

9.525 

1.981 

6.350 

1.194 

OUHUTD 

11.43 

1.321 

7.239 

2.870 

OUHDTD 

11.43 

1.321 

9.271 

0.838 

OUPNTN 

11.43 

1.651 

6.350 

3.429 

OUPUTN 

11.43 

1.651 

8.407 

1.372 

OUHUTU 

11.43 

1.981 

7.239 

2.210 

OUPUTU 

11.43 

1.981 

8.407 

1.041 

OUHDTU 

11.43 

1.981 

9.271 

0.178 

Figure  1.  Top  left:  GCSC  geometric  parameters.  Top  R,ight:  Dimensions  of  tested  inserts.  Bottom  Left:  Film  length 
over  initial  film  height  as  a  function  of  mixture  ratio. ^ 


the  past  decade.  For  example,  Rahman  et  al.  performed  hot  and  cold  experiments  of  a  GCSC  to  obtain  visu¬ 
alizations  of  the  spray  structure  and  measurements  of  the  pressure  drop  across  the  element,^  obtaining  spray 
angle  measurements  but  no  direct  calculations  of  droplet  size.  In  a  preliminary  assessment  of  the  potential 
of  the  GCSC  versus  traditional  coaxial  injectors  Cohn  et  al.  used  a  laser  based  Doppler  interferometer  to 
obtain  measurements  of  droplet  size  and  velocity,  obtaining  a  range  of  diameters  from  3.8  fim  to  440 
during  typical  operating  conditions.®  Soltani  et  al.  investigated  sprays  produced  by  a  liquid-liquid  GCSC 
in  a  non-combusting  environment  to  obtain  measurements  of  average  droplet  size  and  velocity,  noting  that 
both  quantities  achieve  a  self-similar  structure.® 

The  effort  described  below  focuses  on  the  internal  flow  of  a  GCSC  injector,  and  attempts  to  predict  its 
ability  to  perform  over  a  wide  range  of  operating  conditions,  characterizing  the  internal  flow.  In  the  current 
work,  simulations  will  be  used  to  estimate  the  liquid  film  length  and  shape,  particularly  to  predict  whether 
the  film  reaches  the  end  of  the  cup.  Lightfoot  et  al.’s  past  investigations  have  shown  these  film  characteristics 
to  be  very  strongly  correlated  with  the  performance  of  the  injector.^  In  the  current  work,  a  closely-coupled 
two-phase  Eulerian  method  is  investigated  and  used  to  predict  the  models  applicability  to  GCSC  operating 
conditions. 


Experimental  Study 

Experimental  data  was  collected  using  an  injector  body  composed  of  acrylic.  This  acrylic  body  is  modular, 
allowing  differing  injector  geometries  to  be  tested.  Nominally,  the  outlet  diameter  of  the  injector  is  19.05  mm 
(0.75  inches).  The  initial  film  thickness  is  nominally  1.65  mm  (0.065  inches)  with  the  step  separating  the  gas 
and  liquid  having  a  height  of  1.52  mm  (0.06  inches).  Upstream  of  the  modular  acrylic  section  is  a  stainless 
steel  section  which  consists  of  180  mm  of  gas  inlet  with  a  fixed  radius  of  6.35  mm  (0.25  inches).  All  tests 
were  performed  with  atmospheric  back  pressure  using  working  fluids  (simulants)  of  water  and  nitrogen  at 
room  temperature.  The  gas  flow  rates  were  varied  from  0.0187-0.0798  kg/s.  The  liquid  flow  rates  were  varied 
from  0.0236-0.0794  kg/s.  The  momentum  flux  ratio,  defined  using  the  mass  flow  rates  along  with  flow  areas 
based  on  the  initial  film  thickness  for  the  liquid  and  the  average  gas  post  height  -  {rp  +  ro)/2  =  rp-|-(s-|-r)/2 
-  varied  from  around  10  to  around  1100.  Film  length  was  determined  from  video  images.  For  lighting, 
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a  laser  beam  was  split  and  expanded  into  two  sheets  which  were  oriented  180°from  one  another  along  the 
centerline  of  the  injector  or  spray.  A  Vision  Research  Phantom  v7.3  camera  positioned  90°from  the  laser 
sheets  captured  the  images  at  6006  fps  with  an  exposure  time  of  150  microseconds.  A  typical  image,  with 
flow  from  left  to  right,  is  shown  in  Fig.  2.  An  in-house  Matlab  code  was  used  to  determine  the  film  profile 
in  each  image  using  the  change  in  intensity  to  determine  the  film  boundary.  Once  the  film  profile  has  been 
determined,  the  code  uses  input  information  about  the  wall  location  to  determine  film  length.  Because  the 
film  length  is  not  steady,  the  profiles  extracted  from  5000  images  are  averaged  together  to  determine  an 
average  film  length  (and  standard  deviation).  Analysis  shows  that  5000  images  are  sufficient  to  give  good 
statistical  results.  Sample  profiles  are  shown  in  Fig.  3. 


Figure  2.  A  typical  image  from  the  in-cup  video 
is  shown  here.  The  edges  of  the  injector  body  are 
highlighted  including  the  sheltering  lip. 


Figure  3.  Representative  profiles  of  a  long, 
medium  and  short  film.  These  profiles  are  from 
the  geometry  shown  in  Fig.  2  at  momentum  flux 
ratios  of  110,  484  and  823  respectively 


Modeling  Approach 

The  modeling  approach  used  in  the  current  work  recognizes  that  typical  injector  flows  operate  at  ex¬ 
tremely  high  Reynolds  and  Weber  numbers.  The  so-called  E  —  V  model,  originally  proposed  by  Vallet  et 
al.,^  assumes  a  separation  of  the  large  and  small  scales,  allowing  the  bulk  fluid  motion  to  be  directly  simu¬ 
lated  while  the  sub-grid  features  of  the  flow  are  resolved  via  modeling  closures  similar  to  those  used  in  RANS 
turbulence  modeling.  Unlike  Volume  of  Fluid  (VOF),  the  S  —  V  model  does  not  attempt  to  reconstruct  the 
highly  convoluted  interface.  A  VOF  approach  would  be  appropriate  if  one  were  intent  on  direct  numerical 
simulation  and  willing  to  apply  a  micron-level  mesh  to  any  part  of  the  domain  where  small  droplets  might 
exist.  Taking  advantage  of  this  scale  separation  allows  complex  flows  to  be  characterized  with  significantly 
less  computational  resources. 

The  original  model  proposed  by  Vallet  is  founded  on  four  assumptions: 

•  Surface  tension  and  viscosity  act  only  at  small  length  scales,  which  corresponds  to  infinite  Reynolds 
and  Weber  numbers.  This  implies  that  the  large  scale  features  of  the  flow  will  be  dependent  only 
upon  density  variations.  The  surface  tension  and  viscosity  will  only  have  an  effect  upon  the  small  scale 
features. 

•  Although  the  small  scale  velocity  fluctuations  of  the  flow  are  unpredictable  at  the  desired  resolutions, 
the  mean  velocity  field  can  be  predicted  using  standard  turbulence  closures,  such  as  the  classical  two 
equation  k-epsilon  model. 

•  The  dispersion  of  the  liquid  phase  into  the  gaseous  phase  can  be  modeled  via  a  turbulent  diffusion  flux. 

•  The  mean  geometry  of  the  liquid  structures  can  be  characterized  by  tracking  the  mean  surface  area  of 
the  liquid-gas  interface  per  unit  volume 

To  track  the  dispersion  of  the  liquid  phase  an  indicator  function  Y  is  introduced  with  value  1  in  the  liquid 
phase  and  0  in  the  gas  phase.  The  mean  liquid  volume  fraction  of  the  fluid  is  then  given  as  Y,  and  the  mean 
mass  averaged  fraction  is  defined  as  Y  =  ^.  The  transport  equation  for  Y  then  takes  the  form: 

dpY  dpUiY  dpu'^Y' 

where  u'  denotes  the  turbulent  fluctuations  in  velocity  and  Y'  denotes  turbulent  fluctuations  in  volume 
fraction. 

The  value  of  Y  is  related  to  p  by 
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Treating  the  transport  equation  in  this  way  allows  the  momentum  equation  to  be  solved  using  the  single 
equation: 


dpuj  dpuiUj  dp  dpuiUj  ,  , 

dt  dxi  dxj  dxi 

The  turbulent  diffusion  liquid  flux  u'^Y'  captures  the  effect  of  the  relative  velocity  between  the  two  phases. 
Models  for  both  the  turbulent  diffusion  liquid  flux  in  the  transport  equation  and  the  turbulent  diffusion  flux 
in  the  momentum  equation  must  be  formed  to  obtain  closure.  The  further  assumption  that  both  the  gas 
and  liquid  density  is  constant  closes  the  system.  Most  previous  implementations  of  the  YY  model  use  a 
classical  compressible  k  —  e  model  to  close  the  turbulent  diffusion  term,  although  Demoulin  et  al.  proposed 
a  modification  of  the  transport  equation  to  capture  the  effects  of  a  large  density  gradient  on  the  production 
of  kinetic  energy.® 

With  the  large  scale  flow  features  described,  a  transport  equation  for  the  evolution  of  the  interfacial  surface 
area  determines  the  small  scale  behavior  of  the  flow.  There  have  been  several  different  forms  suggested  for 
this  equation.  A  rigorous  derivation  of  the  equation  from  first  principles  has  been  performed  by  Morel.®  The 
transport  equation  of  S  proposed  by  Vallet^®  takes  the  same  form  as  Morel’s  generalized  derivation,  having 
terms  capturing  the  effects  of  convection,  diffusion,  surface  stretching  and  coalescence. 


?  +  ^  +  — S  -  — E® 

dt  dXj  dXj  \  dXj  J  Tprod  Tdestr 


(8) 


Tprod  is  a  time  scale  representing  the  rate  at  which  turbulence  creates  surface  area,  and  Tdestr  is  a  velocity 
scale  corresponding  to  the  rate  at  which  collision  and  coallesence  destroy  surface  area.  A  detailed  derivation 
of  these  terms  can  be  found  in  Vallet’s  original  description  of  the  model. ^  Since  the  model  was  initially 
published,  several  alternative  forms  of  this  equation  have  been  proposed  to  better  resolve  differences  in 
primary  and  secondary  atomization.^^  Knowledge  of  E  allows  the  prediction  of  the  Sauter  mean  diameter 
and  the  droplet  number  density  via 


^32  = 


6pY 

PiY 


(9) 


and 


367r;o2f2 


(10) 


A  transport  equation  having  this  form  predicts  an  equilibrium  value  of  E  when  Ygq  = 


'^destr 

'^prod 


Diffusion  liquid  flux  closure 

The  ability  of  the  model  to  effectively  predict  the  structure  of  the  liquid  film  is  governed  by  the  closure 
of  the  diffusion  flux  in  the  volume  fraction  transport  equation.  Vallet  et  al.  derived  an  algebraic  relation 
for  this  flux  by  first  considering  the  following  transport  equation  for  the  mass  fluctuation  Y'  derived  via  the 
Reynolds  averaging  process.^® 


dpu[Y'  dpUju[Y' 
dt  dxj 


dui 


dxi 


y./  dp  ai  ~ 

Y  - - p—u\Y' 

dXi  Tt 

(11) 


The  first  term  on  the  right  hand  side  corresponds  to  the  diffusion  of  mass  fluctuations,  the  following  three 
terms  to  production,  and  the  final  term  to  destruction. 

Equation  11  assumes  that  can  be  expressed  as 


,dY'  _ai  '—XT,  V  ^YrYdui  lYxri  dp 

P  =  -p—u'^Y'  +  j^pu'jY'- - h7  E  — 

OXi  Tt  OXj  OXi 


(12) 
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where  Tt  is  the  integral  turbulence  timescale  and  Ui,  7^  and  7'^  are  constants.  Assuming  dominance  by 


the  production  and  destruction  terms  yields  the  following  algebraic  closure. 


pu\Y'  =  -  — 

ai 


where  Y'  is  defined  exactly  as 


dui 


JY\  dp 


dxj 


(13) 


(14) 


1  1 
.Pi  Pg, 

Under  the  further  assumption  that  diffusion  is  dominated  by  the  first  term  on  the  right  hand  side  of 


Equation  13,  this  results  in  the  classical  Pick’s  law  of  diffusion 

pt  dY 


(15) 


Demoulin  proposed  the  following  term  to  capture  the  enhancement  of  Rayleigh-Taylor  instability  growth 
caused  by  the  presence  of  a  large  density  gradient. 


pu'iY'  =  p 


f )y(i-y) 


Sc 


Pg  Pi 


dY 

dxi 


(16) 


The  inclusion  of  this  term  has  proven  successful  in  resolving  the  primary  liquid  core  size  for  both  diesel^^ 
and  coaxial®  injections. 


Implementation 

The  Yi  —  Y  model  was  implemented  using  the  OpenFOAM  software  framework. The  implementation 
uses  a  fully  parallelized  finite-volume  method  generalized  to  arbitrary  convex  polyhedral  cells  in  two  and 
three  dimensions.  Discretization  methods  are  selected  at  run-time  and  include  over  thirty  choices  of  flux 
limiting. 

A  segregated  pressure-based  PISO  method^®  has  been  implemented  to  update  variables  at  each  timestep 
by: 

•  Solving  the  continuity  equation  using  the  volumetric  flux  from  the  previous  timestep 

•  Solving  the  species  transport  equation  using  the  mass  flux  predicted  by  the  continuity  equation 

•  Solving  the  momentum  equation  to  predict  the  current  volumetric  flux 

•  Iteratively  solving  the  pressure  equation  until  the  velocity  field  satisfies  the  continuity  equation 

•  Updating  continuity  and  species  transport  equations 

•  Either  iterate  to  correct  for  coupling  between  species  transport  and  pressure  equation,  or  procede  to 
next  timestep 

To  derive  the  pressure  equation,  the  single  phase  incompressible  pressure  equation  typically  implemented 
in  OpenFOAM®"*  was  altered  to  account  for  variable  density  effects.  The  semi-discretized  form  of  the  mo¬ 
mentum  equation  is  defined  as 


Op  Up  =  H(U)  -  Vp  (17) 

where  H(u)  consists  of  two  parts:  the  advection  terms  which  includes  the  matrix  coefficients  of  neigh¬ 
boring  cells  multiplied  by  the  corresponding  velocities  and  the  source  term,  accounting  for  all  point  sources 
in  the  momentum  equations  (such  as  gravity,  body  forces,  etc.). 

H(U)  =  ^a^U^  +  ^  (18) 

N 
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Equation  17  is  rearranged  and  interpolated  to  faces.  Taking  the  divergence  of  this  gives 


V  .  (U,),  =  V  .  (1h(U),)  -  V  .  (Ivp) 


(19) 


The  effect  of  variable  density  is  then  coupled  to  the  pressure  equation  via  Equation  6,  the  continuity 
equation  and  the  chain  rule 


Solving  the  equation 


V  •  (Up)/ 


1  dp  DY 
pdY  Dt 


pu\Y' 


V  • 


pu'iY' 


and  updating  the  fluxes  will  give  a  velocity  field  that  satishes  the  continuity  equation. 


(20) 


(21) 


Results 

The  conditions  and  geometry  for  this  particular  study  are  discussed  in  Lightfoot  et  al.^  Specifically,  the 
ONPNTN  geometry  was  considered  under  the  operating  conditions  denoted  in  Table  1.  The  computational 
boundary  conditions  were  set  by  matching  liquid  flow  rate  and  gas  flow  rate,  assuming  a  uniform  incoming 
velocity  profile,  and  a  fixed  total  pressure  condition  at  the  injector  outlet.  It  is  important  to  note  that  for 
several  of  the  higher  momentum  ratio  cases,  the  gas  flowrates  yield  flow  velocities  in  the  transonic  regime, 
rendering  the  assumption  of  incompressibility  questionable.  Experimental  measurements  of  the  total  pressure 
at  the  inlet  suggest  that  choking  is  not  occuring,  but  the  difference  in  gas  density  could  substantially  alter 
the  momentum  ratio  at  the  interface.  To  obtain  an  estimate  of  the  effect  of  compressibility,  a  separate 
single-phase,  turbulent,  compressible  flow  solver  was  used  to  compare  the  predicted  velocity  and  density 
fields.  This  solver  obtained  a  range  of  densities  from  0.98  —  1-4^,  suggesting  a  15%  possible  difference  in 
gas  momentum  at  the  interface. 


Table  1.  Simulation  conditions 


Case  name 

Pgas 

Pliq 

Ugas^in 

Uliq^in^z 

Uliq^in^O 

Pout 

kq 

kq 

m 

s 

m 

s 

m 

s 

kPa 

ONPNTN-al 

1.2 

999 

300 

0.35 

4.2 

101.325 

ONPNTN-bl 

1.2 

999 

153 

0.52 

6 

101.325 

In  the  computation,  the  domain  is  in  two-dimensional  polar  coordinates  (r,  2),  with  r  >  0,  with  a  two- 
dimensional  mesh  containing  eleven  thousand  cells.  The  assumption  of  axisymmetric  flow  was  validated 
by  performing  a  3D  simulation  of  a  90°  sector  of  the  full  geometry.  While  the  experimental  photographs 
did  exhibit  asymmetry  between  the  top  and  bottom  film  prohles,  this  has  been  attributed  to  experimental 
conditions  and  is  discussed  by  Schumaker  et  al.^® 

Turbulence  was  modeled  using  both  a  classical  two  equation  fc  —  e  model  and  the  realizable  k  —  e  model 
developed  by  Shih.^®  All  previous  implementations^dii  17,  is  Y  —  Y  model  have  used  the  k  —  e  model 

because  of  its  stability  and  ability  to  resolve  mixing  layers,  which  was  well  suited  for  the  plain  orifice  and 
coaxial  arrangements  previously  considered.  The  k  —  e  model  is  well  known  to  be  unable  to  accurately 
predict  the  recirculation  zone  following  a  backward  facing  step.  While  any  turbulence  model  based  on  the 
Boussinesque  assumption  will  have  difficulty  capturing  the  effect  of  swirl,  the  realizable  k  —  e  model  should 
be  capable  of  more  accurately  resolving  the  strength  of  the  recirculation  zone,  which  is  assumed  to  be  the 
dominant  turbulent  mechanism  affecting  the  film  profile. 

Simulations  were  run  using  the  standard  k  —  e  model  for  the  ONPNTN-bl  case  using  the  default  value 
for  the  Schmidt  number  of  0.9.  Figure  4  shows  the  resulting  pressure  and  volume  fraction  profiles.  Figure  5 
shows  an  instantaneous  photograph  of  the  injector  under  the  same  operating  conditions.  The  experimental 
photograph  suggests  that  the  liquid  should  be  pushed  up  against  the  backward  facing  step  by  the  recirculation 
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zone,  and  the  lack  of  this  feature  suggests  that  the  model  is  underpredicting  the  turbulent  dispersion.  Also 
visible  in  Figure  4  is  a  static  pressure  gradient  caused  by  the  centripetal  acceleration  present  in  the  swirling 
liquid  him. 


p(Pa) 

101325. 

1.  106644. 

Volume  Fraction 

0.0  0.^^^^^^  0.75  1.0 

Figure  4.  Static  pressure  and  volume  fraction  profile  for  ONPNTN-bl  case 


Preliminary  results  using  Demoulin  et  al’s  modihed  closure  (Equation  16)  further  underestimated  the  him 
length.  This  is  attributed  to  the  fact  that,  for  the  swirling  liquid  him  under  consideration,  the  centripetal 
acceleration  acting  upon  the  interface  has  a  stabilizing  effect  that  supresses  the  growth  of  the  Rayleigh-Taylor 
type  instability. 

As  a  hrst-order  approximation  of  the  reduced  turbulent  dispersion  caused  by  the  swirl,  the  Schmidt 
number  was  increased  to  determine  its  ehect  upon  the  liquid  him  prohle.  Figure  6  demonstrates  that 
increasing  the  Schmidt  number  from  the  standard  value  of  0.9  up  to  a  value  of  20  allows  the  model  to 
accurately  reproduce  the  experimental  data. 

Figure  7  gives  a  similar  parametric  study  of  the  Schmidt  number  when  the  realizable  k  —  e  turbulence 
closure  is  used  instead.  For  this  case,  a  Schmidt  number  of  only  15  provides  a  more  accurate  representation 
of  the  him  prohle.  This  suggests  that  the  standard  k  —  e  model  is  overpredicting  the  liquid  dispersion,  as 
expected. 

To  characterize  the  generality  of  these  results,  a  case  with  a  signihcantly  shorter  him  length  was  inves¬ 
tigated.  Figure  8  shows  an  experimental  photograph  of  the  him  prohle.  For  this  case,  the  light  scattering 
at  the  interface  is  substantially  more  diffuse,  suggesting  the  presence  of  more  turbulent  mixing.  Figures  9 
and  10  compare  the  predicted  liquid  prohles  to  the  time-averaged  experimental  data.  Both  simulations  show 
good  initial  agreement,  but  are  abruptly  cut  short  as  a  recirculation  zone  forms  at  the  end  of  the  liquid 
him  (Figure  11).  As  expected,  the  realizable  k  —  e  model  comes  closer  to  predicting  the  experimental  result 
but  also  fails  to  accurately  resolve  the  correct  behavior  near  the  recirculation  zone.  In  Beheshti  et  al.’s 
assessment  of  the  E  —  F  model,  they  noted  that  in  its  current  form  the  model  is  incapable  of  accurately 
resolving  hows  with  areas  of  large  recirculation.  In  this  scenario,  the  fundamental  assumption  that  the  how 
is  acting  at  large  Reynold’s  and  Weber  number  is  violated,  as  the  velocity  within  the  recirculation  zone  is 
substantially  less  than  the  free-stream  velocity. 

Conclusions 

Resolving  the  complex  huid  mechanics  occuring  within  the  cup  of  a  GCSC  presents  numerous  experi¬ 
mental  and  modeling  challenges.  The  correspondence  of  both  the  experimental  and  computational  results 
is  consistent  with  the  him  extent  being  determined  by  turbulent  transport.  For  CFD  simulation,  rather 
than  trying  to  resolve  the  wide  variety  of  scales  present  in  the  cup,  the  resolved  scales  can  be  treated  with 
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Figure  5.  Experimental  photograph  of  liquid  profile  for  ONPNTN-bl  case 


Figure  6.  Iso-contours  of  Y  —  0.5  compared  to  experimental  profiles  for  standard  k  —  e  closure 
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Figure  7.  Iso-contours  of  =  0.5  compared  to  experimental  profiles  for  realizable  k  —  e  closure 


Figure  8.  Experimental  photograph  of  liquid  profile  for  ONPNTN-al  case 
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Figure  9.  Iso-contours  of  y  =  0.5  compared  to  experimental  profiles  for  standard  k  —  e  closure 


X  10'^ 


Figure  10.  Iso-contours  of  y  =  0.5  compared  to  experimental  profiles  for  realizable  k  —  e  closure 
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Figure  11.  Instantaneous  film  profile  and  streamlines  show  the  presence  of  a  recirculation  zone  for  the  ONPNTN-al 
case 


conservation  laws  and  the  unresolvable  turbulence/interface  interaction  can  be  represented  by  models.  While 
the  model  was  able  to  predict  low  momentum  ratio  film  profiles,  it  failed  to  predict  film  structure  when  the 
gas  flowrate  was  increased  and  recirculation  zones  became  more  pronounced.  This  shortcoming  has  not  yet 
been  observed  in  other  injector  arrangements.  For  this  model  to  be  used  as  a  reliable  predictor  of  film  length, 
higher  order  modeling  closures  for  the  turbulent  flux  must  be  used  to  resolve  the  complex  recirculation  zone 
interactions. 
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